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HIGHLIGHTS 


•  Generalized  moving  boundary  model  for  lithium  ion  cells  with  phase  change  electrode. 

•  Numerical  implementation  is  obtained  using  a  pseudo  2D  approach. 

•  Model  captures  tangential  phase  front  movement  and  phase  separation  in  LiFeP04. 

•  Simulation  shows  charge-discharge  asymmetry  and  path  dependency  in  LiFeP04. 
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Lithium  ion  cells  with  electrode  materials  that  undergo  phase  transitions,  like  LiFeP04,  have  unique 
charge— discharge  characteristics.  In  this  work  a  generalized  framework  of  moving  boundary  approach  is 
proposed  to  model  the  path  dependent  charge-discharge  response  of  porous  electrodes  that  exhibit 
multi-phase  coexistence.  Using  the  Landau  transformation  the  governing  equation  in  moving  coordinate 
is  transformed  to  fixed  coordinate  and  a  suitable  path  dependent  algorithm  is  devised  and  is  imple¬ 
mented  in  a  multi-physics  environment.  Simulation  results  show  that  tangential  propagation  of  the 
phase  front,  often  seen  in  experiments,  can  be  addressed  by  this  model.  Incorporation  of  multi-phase 
diffusion  predicts  the  characteristic  phase  separation  in  LiFeP04  particles.  The  proposed  model  suc¬ 
cessfully  captures  the  charge— discharge  asymmetry  of  LiFeP04  based  cells.  Efficacy  of  the  proposed 
approach  to  model  the  path  dependence  of  cells  with  phase  change  electrodes  is  demonstrated  by 
simulating  the  response  of  LiFeP04/graphite  cell  subjected  to  a  charge-discharge  pulse.  Numerical 
studies  are  performed  to  study  the  effect  of  important  model  parameters  to  enhance  cell  design  and  to 
bring  out  unique  features  in  the  cell  response  due  to  multi-phase  coexistence. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Lithium  iron  phosphate  (LiFePCU,  LFP)  is  increasingly  deployed 
as  positive  electrode  material  in  lithium  ion  cells  owing  to  its  high 
thermal  stability,  cost  effectiveness,  non-toxic  nature  and  long  cy¬ 
cle  life  [1—3].  Given  these  advantages,  it  however  suffers  from 
limitations  [4-6]  due  to  low  electronic  conductivity  resulting  in 
higher  Ohmic  drop  in  the  electrode  and  low  rate  capability  leading 
to  reduced  active  material  utilization.  Several  research  efforts  over 
the  last  decade  have  enabled  significant  improvements  by 
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enhancing  properties  of  the  material  through  suitable  modifica¬ 
tions.  Some  of  the  prominent  ideas  pertain  to  coating  of  LFP  par¬ 
ticles  with  carbon  to  improve  electronic  conductivity  [7,8],  decrease 
in  particle  size  for  reduction  in  transport  losses  [9,10],  use  of  LFP 
nano-crystals  to  increase  the  surface  area  [11],  solid-solution 
doping  by  supervalent  metals  [12]  etc.  Because  of  these  de¬ 
velopments,  LFP  based  cells  have  found  applications  in  mobile, 
automotive  and  other  industries. 

The  wide  spectrum  of  applications  of  LFP  based  cells  encourages 
electrochemical  modeling  to  understand  the  underlying  cell  pro¬ 
cesses.  In  addition  to  optimized  cell  design,  electrochemical  models 
enable  robust  battery  management  systems  (BMS)  development 
for  on-board  applications.  A  critical  challenge  in  developing  BMS 
algorithm  for  cells  with  LFP  based  electrodes  is  the  non-monotonic 
dependence  of  open  circuit  potential  on  capacity.  This  behavior  is 
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also  observed  in  other  electrode  materials,  mainly  due  to  the  un¬ 
derlying  crystallographic  solid-solid  phase  transformation. 
Therefore,  modeling  phase  change  is  an  active  area  of  research  that 
enables  development  of  a  complete  electrochemical  model  [8]  for 
lithium  based  cells.  These  models  are  proposed  at  various  length 
scales  varying  from  crystalline  level  to  continuum  level.  Existing 
approaches  in  electrochemical  modeling  of  LFP  electrodes  can  be 
classified  into  three  categories  viz.,  (a)  shrinking  core  models  [13- 
18],  (b)  phase  field  models  [19-23],  and  (c)  resistive  reactant 
models  [24-28].  A  brief  overview  of  these  models  is  provided  in 
Table  1. 

An  electrochemical  model  based  on  a  core-shell  approach  was 
first  proposed  for  LFP  electrode  by  Srinivasan  and  Newman  [14]  and 
later  developed  [15]  for  an  LFP/C  cell.  As  mentioned  in  Table  1,  this 
approach  captures  the  discharge  behavior  by  considering  interca¬ 
lation  coupled  mass  transport  due  to  diffusion  in  a  Li  rich  p  phase 
and  movement  of  two-phase  interface.  In  this  framework  however, 
the  diffusion  in  Li-poor  phase  (a)  is  neglected.  The  underlying 
construct  of  this  model  assumes  diffusion  in  a  particle  causing 
lithium  flux  movement  normal  to  the  phase  boundary.  On  the 
contrary,  experimental  studies  [1]  reveal  Li  movement  in  a 
tangential  direction.  The  diffusion  under  phase  transition  is 
modeled  using  sharp  interface  model  which  predicts  a  fast  Li 
transfer  across  the  phase  boundary  between  a  and  (3  phase.  Fol¬ 
lowed  by  this  work,  a  more  sophisticated  version  of  the  core— shell 
model  in  which  diffusion  in  both  phases  of  LiCo02  cathode  [13]  is 
developed.  Also,  Wang  and  co-workers  in  Ref.  [16],  proposed  a 
model  for  LFP,  where  the  interface  movement  between  the  two 
coexisting  phases  was  considered  to  be  sluggish  by  incorporating  a 
mobility  term  in  the  interface  movement  equation.  Kasavajjula  [17] 
further  enhanced  this  model  by  considering  the  Li  diffusion  in  both 
phases  together  along  with  the  interface  mobility.  Recently,  a 
multi-scale  approach  proposed  by  Dargaville  and  Farrell  [18] 
employed  the  shrinking-core  model  of  Srinivasan  and  Newman 
[14]  to  represent  the  three  size  scales  present  in  a  porous  LFP 
cathode,  namely  crystals,  agglomerates,  and  the  porous  cathode 
itself. 

Above  developments  show  that  the  core-shell  approach  pro¬ 
vides  a  mathematically  consistent  framework  to  incorporate  effects 
of  phase  transitions  in  a  complete  electrochemical  model.  The 
existing  models  [13-18]  do  not  however  account  for  factors  such  as 
diffusion  in  the  Li-poor  LFP  phase  (a),  Ohmic  drop  due  to  carbon 
coating  to  model  inhomogeneous  current  distribution  etc.  As  a 
result,  the  application  of  the  shrinking  core  approach  is  limited  to 
represent  discharge  behavior  only.  Safari  and  Delacourt  [26,27] 
recently  proposed  a  reduced  order  electrochemical  model  using 
the  resistive  reactant  concept  without  incorporating  the  details  of 
phase  transformation  in  LFP  electrode.  It  successfully  captures  the 
charge— discharge  asymmetry  but  does  not  represent  the  features 
due  to  phase  transformation  in  terms  of  tangential  front  propaga¬ 
tion,  phase  separation  and  path  dependence.  Hence,  development 
of  a  suitable  electrochemical  approach  that  models  the  multi-phase 
coexistence  during  charge  and  discharge  consistently  is  highly 
desirable.  Considering  the  simple  and  robust  framework  of 
shrinking  core  model  an  extension  of  this  approach  to  represent 
multiple  charge-discharge  cycle  is  yet  to  be  attempted.  Such  a 
general  framework  would  also  enable  detailed  investigation  of 
charge— discharge  asymmetry  and  path  dependence  [29] 
frequently  reported  for  LFP  electrodes. 

In  the  present  work  a  general  electrochemical  framework  is 
proposed  that  identifies  the  mechanisms  governing  voltage 
response  in  cells  with  phase  change  electrodes.  In  the  following 
sections,  a  generalized  path  dependent  moving  boundary  model  is 
developed  and  subsequently  modified  as  a  cell  design  tool.  Section 
2  describes  the  phase  transformation  processes  in  LFP  which  brings 
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VII  VI 


a  -  Li  poor  phase 
p  -  Li  rich  phase 


Fig.  1.  Schematic  diagram  describing  the  various  stages  of  lithium  intercalation  and  de-intercalation  in  LFP  during  charging  and  discharging  processes  based  on  core-shell 
construct. 


out  the  physics  of  path  dependence  under  charge-discharge  cycles. 
Details  of  the  governing  equations  of  the  proposed  framework  are 
provided  in  Section  3.  In  Section  4  the  path  dependent  algorithm  is 
devised  and  steps  involved  in  numerical  implementation  are  briefly 
discussed.  Section  5  discusses  the  results  obtained  from  the  nu¬ 
merical  simulations  of  proposed  model.  This  section  starts  with  the 
proposed  approach  being  validated  with  results  [15]  reported  in 
literature.  Subsequently  the  charge-discharge  asymmetry  and 
path  dependency  in  LFP/C  cells  are  investigated  using  several  case 
studies.  Key  insights  from  the  simulations  that  can  be  of  relevance 
for  optimal  cell  design  are  also  proposed. 

2.  Phase  transformation  processes  in  LFP  particle 

LFP  is  naturally  available  as  a  phosphate  mineral  and  is  an 
olivine-type  material.  It  has  triphylite  structure  belonging  to 
orthorhombic  system  of  Bravais  lattices.  Ability  to  intercalate  and 
de-intercalate  lithium  ion  at  desirable  range  of  voltage  and  current 
enables  application  of  this  material  as  positive  electrode  material  of 
secondary  lithium  ion  cells.  Interestingly  LiFePCU  is  also  known  to 
undergo  crystallographic  structural  phase  transformation  from 
lithiated  triphylite  structure  (denoted  by  (3)  to  de-lithiated  heter¬ 
osite  (denoted  by  a)  structure  during  reversible  extraction  of 
lithium.  Padhi  and  co-workers  [1]  discuss  the  details  of  the  two- 
phase  coexistence  and  crystallographic  structure.  Further  Dodd 
[30]  and  Ramanna  [31]  use  advanced  material  characterization 
techniques  to  provide  useful  insights.  As  discussed  by  Padhi  [1  ],  the 
open  circuit  voltage  curve  of  Li^FePCU  is  constant  over  a  large  range 
of  values  of  the  intercalate  concentration  x.  It  can  be  ascertained 
using  the  Gibb’s  phase  rule  that  two  phase  coexistence  drives  the 
underlying  intercalation/de-intercalation  process  resulting  in  a 
non-monotonic  open  circuit  voltage.  An  outcome  of  the  multi¬ 
phase  coexistence  in  this  material  is  a  path  dependent  response. 


One  of  the  approaches  to  represent  this  multi-phase  coexistence  is 
a  sharp  interface  (or  the  core-shell)  model. 

A  schematic  of  the  phase  transformation  process  in  LFP  is 
shown  in  Fig.  1  [29].  Various  stages  (I-X)  during  lithiation/deli- 
thiation  in  a  two-phase  material  of  spherical  geometry  is  illus¬ 
trated  in  Fig.  1  using  the  core-shell  approach.  During  discharge 
process,  as  current  is  drawn  from  the  cell,  lithium  intercalates  into 
the  de-lithiated  LFP  particles  that  constitute  the  cathode,  as 
shown  in  stage  I  and  stage  II.  Onset  of  phase  transformation  from 
de-lithiated  a  to  lithiated  [3  begins  after  reaching  a  critical  level  of 
concentration  (c%q).  This  leads  to  phase  separation  with  a  shell  of 
lithium  rich  [3  phase  forming  over  the  core  of  lithium  poor  a  phase, 
as  shown  in  stage  II  ->  III  ->  IV.  Stage  III  and  stage  IV  of  schematic 
Fig.  1  show  that  during  the  course  of  phase  transformation  the 
phase  separation  occurs  and  core  shrinks  due  to  gradual  move¬ 
ment  of  the  two  phase  interface.  In  phase  change  materials  the 
mass  transfer  is  directly  related  to  movement  of  interface  which  is 
in  turn  is  related  to  constant  discharge  voltage.  Phase  trans¬ 
formation  ceases  as  the  interface  reaches  the  center  (Stage  V  in 
Fig.  1)  phase  transformation  stops  and  the  material  behaves  as 
single  phase  system.  Diffusion  in  pure  (3  phase  continues  from 
stage  V  to  stage  VI  and  the  charge  transfer  reaction  stops  on 
successive  completion  of  filling  of  all  de-lithiated  sites,  and  this 


Fig.  2.  Schematic  showing  formation  of  multiple  cores  and  shells  in  the  event  of 
charging  of  incompletely  discharged  particle. 
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marks  end  of  discharge.  When  the  completely  discharged  or 
lithiated  particle  is  charged  (stage  VI),  the  above  described  pro¬ 
cess  repeats  (stage  VI  to  stage  X).  The  difference  during  charging  is 
that  a  phase  forms  the  shell  and  [3  phase  forms  the  core.  In 
contrast  to  discharging,  the  single  ((3)  phase  exists  for  a  larger  SOC 
range.  During  the  SOCs  when  both  the  phases  coexist,  the  charge 
transfer  reaction  occurs  in  lithium  deficient  a  phase.  From  this 
discussion,  it  is  clear  that  the  diffusion  and  charge  transfer  process 
predominantly  occur  in  different  materials  during  charging  and 
discharging  [29]  and  hence  the  response  different  between  the 
processes.  This  qualitatively  explains  the  existence  of  asymmetric 
behavior  in  LFP  particle  which  is  in  general  a  feature  [29,32]  of  any 
phase  changing  electrode.  More  importantly  the  phase  trans¬ 
forming  material  also  shows  path  dependent  response  due  to 
multiple  phase  coexistence.  This  is  further  explained  using  a 
construct  wherein  the  lithiation  process  during  discharge  is  not 
complete  and  the  charging  starts.  Under  such  situation  shell  of 
lithium  poor  phase  forms  (Fig.  2)  which  can  lead  to  multiple  core 
shell  formation  and  eventually  path  dependent  response.  Multiple 
core  shell  formation  will  also  accentuate  the  phase  separation  and 
phase  inhomogeneity.  Quantitative  explanation  of  this  behavior  is 
attempted  in  the  present  work  by  devising  and  numerically 
solving  a  generalized  moving  boundary  model  in  the  pseudo  2-D 
framework  of  porous  electrode  theory.  Details  of  the  generalized 
governing  equations  are  given  in  the  next  section. 


V-(<7SV0S)  -  asin  =  0 


(3) 


It  may  be  noted  that  the  electrodes  are  treated  as  homogenized 
continuum  with  effective  properties  given  by  Bruggeman’s  rela¬ 
tionship.  The  model  equations  are  adopted  for  three  different  zones 
in  the  cell  viz.,  anode,  separator  and  cathode  based  on  which  the 
boundary  conditions  are  prescribed.  The  quantity  cs  appearing  in 
Eq.  (5)  is  the  lithium  concentration  at  the  surface  of  the  particle, 
obtained  by  solving  solid  phase  mass  transport  equation  given  by 
Fields  law.  The  mass  transport  equation  needs  to  be  modified  for 
materials  undergoing  phase  transformation,  as  in  the  present  case 
of  LFP  particles.  Using  the  sharp  interface  model  the  interface 
coupled  diffusion  equation  [38]  in  spherical  coordinates  for  mate¬ 
rials  undergoing  phase  transition  can  be  expressed  as: 
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3.  Electrochemical  model  for  charge-discharge  with  phase 
change  electrode 

Continuum  description  of  coupled  electrochemical  mass  and 
charge  transport  in  lithium  ion  cells  is  obtained  using  porous 
electrode  theory.  Newman  [33,34]  pioneered  this  approach  and  it 
has  led  to  development  of  an  efficient  mathematical  tool  that 
captures  lithium  diffusion  dynamics  and  charge  transfer  kinetics  to 
predict  current/voltage  response  of  a  battery.  This  approach  is  used 
to  provide  design  guidelines  using  theoretical  underpinning  of 
thermodynamics,  kinetics,  and  transport  across  electrodes. 

In  the  Refs.  [33-37]  details  of  the  model  based  on  porous 
electrode  theory  can  be  found.  At  the  outset,  the  important  as¬ 
sumptions  made  in  porous  electrode  theory  are  listed  below: 

1.  Electrode  is  represented  as  a  superposition  of  active  material, 
filler,  and  electrolyte  that  coexist  at  every  point. 

2.  Particles  of  the  active  material  are  modeled  as  spheres. 

3.  The  electrode  phase  is  coupled  to  the  electrolyte  phase  by  the 
charge  transfer  reaction. 

4.  All  phases  are  considered  to  be  electrically  neutral. 

Based  on  these  assumptions  the  electrochemical  model  for 
lithium  ion  cell  is  mathematically  described  by  the  following  gov¬ 
erning  partial  differential  equations  (PDEs).  These  consist  of  mass 
transport  equation  (Eq.  (1 )),  modified  Ohm’s  law  for  electrolyte  (Eq. 
(2)),  Ohm’s  law  for  the  electrode  (Eq.  (3))  and  Butler-Volmer 
equation  (Eqs.  (4)  and  (5))  which  describes  reaction-diffusion  ki¬ 
netics  that  couple  the  concentration  and  potential  fields  of  elec¬ 
trode  and  the  electrolyte. 
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C  =  C0|t=0  V(x,r) 

where  rj"  ~ p  denotes  the  interface  position,  us_>p  denotes  the  ve- 
locity  of  the  interface  between  source  (s)  and  product  (p)  phases 
that  undergo  phase  transition.  The  step  function  5,  is  defined  to  be 
unity  at  r^p.  At  the  beginning  of  discharge  of  LFP  cathode,  the 
particle  gets  lithiated  and  hence  the  transformation  is  from  a  ->  (3. 
In  other  words  a  is  the  source  and  (3  the  product.  Consequently, 
subscript  T  in  the  diffusivity  coefficient  can  either  be  ‘s’  or  ‘p’ 
depending  upon  the  regime  of  phase  coexistence.  General  expres¬ 
sion  for  interface  boundary  velocity  is  obtained  by  the  Stefan 
condition  based  on  the  mass  balance  at  the  interface  and  is  given 
as: 


,s^p 
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The  sharp  interface  model  imposes  the  following  Dirichlet 
boundary  condition  at  the  interface 


c  =  cpeq  at  r  =  rpp(t) 


(10) 


-  V  •  (Df f Vce)  -  L-2+  asin  =  0  ( 1 ) 

+  asin  =  0  (2) 


The  criterion  for  phase  transformation  is  given  as  c  >  cseq.  It  is  to 
be  noted  that  the  boundary  condition  (Eq.  (10))  that  separates  the 
two  phase  is  not  fixed  in  space.  The  interface  variable  r2s_*p  (t), 
which  denotes  position  of  the  interface,  evolves  with  time  as  given 
by  Eq.  (9).  The  set  of  equations  required  to  account  for  phase 
transitions  in  electrode,  cathode  in  the  present  case,  are  given  by 
Eqs.  (6)— (10).  The  nature  of  equations  is  kept  general  enough  to 
account  for  forward  and  reverse  transformation  during  charge/ 


(c  >  ceq)  Charging  (c  <  c?q) 
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S  II  II 


Table  3 

Geometric  and  material  properties  for  LFP/C  cell  used  in  present  simulations. 


Parameter 

Symbol 

Natural 

graphite 

Iron-phosphate 

Electrode  thickness  [15]  (pm) 

5 

50 

75 

Electrode  porosity  [15] 

£ 

0.33 

0.27 

Average  particle  radius  [15] 

R 

11  pm 

52  nm 

Density  [15]  (kg  m-3) 

P 

3370 

3600 

Diffusion  coefficient  [15] 

Ds 

9  x  10~14 

8  x  1018 

(solid  phase)  m2  s  1 

Reaction  rate  constant  [15] 

k 

3  x  10"5 

3  x  1017 

Matrix  phase 

(Ts 

100 

5  x  10"3 

conductivity  [15] 

(Sm-1) 

Maximum  concentration 

-max 

cs 

27,760 

20,950 

[15]  (mol  m  3) 

Parameter 

Separator  thickness  [15]  (pm) 

Is 

Value 

25 

Separator  porosity  [15] 

£sep 

0.55 

Separator  density  [15]  (kg  m-3) 

P  sep 

900 

Electrolyte  density  [15]  (kg  m-3) 

ki 

1200 

Equilibrium  concentration 

Ceq 

0.048cmax 

for  a  [16] 

Equilibrium  concentration 

cL 

0.8920cmax 

for  3  [16] 

discharge  or  simultaneous  charge-discharge.  The  specific  details  of 
these  equations  (Eqs.  (6)— (10))  and  the  coupling  with  Eqs.  (1)— (5) 
using  the  pseudo  2-D  approach  is  discussed  in  the  next  section.  For 
electrodes  without  phase  transitions,  however  (anode  in  the  pre¬ 
sent  case)  Eq.  (6)  without  the  interface  movement  term  along  with 
boundary  and  initial  condition  as  in  Eqs.  (7)  and  (8)  respectively 
characterizes  the  diffusion  process. 

4.  Path  dependent  model  for  LFP/C  cell  and  numerical 
implementation 

Interface  coupled  diffusion  problem  as  described  in  the  earlier 
section  (Eqs.  (6)— (10))  falls  under  the  general  class  of  inverse 
mathematical  problem  called  Stefan  [39]  problem.  These  kinds  of 
problems  arise  in  several  phase  transition  phenomenon  like  so¬ 
lidification  of  molten  metal  [40],  crystal  growth  [38]  and 
martensitic  transformation  [41  ]  during  steel  forming.  In  order  to 
numerically  solve  Eqs.  (6)— (10)  several  methods  and  algorithms 
[13,39]  are  proposed.  These  methods  can  be  classified  as  fixed  mesh 
method  and  moving  mesh  method.  Moving  mesh  method  [42]  uses 
ALE  (Arbitrary  Lagrangian-Eulerian)  approach  and  allows  for  mo¬ 
tion  of  the  boundary.  The  mesh  motion  is  defined  in  Eulerian  co¬ 
ordinates  and  is  independent  of  motion  of  the  material,  defined  in 
Lagrangian  coordinates.  This  method,  however,  requires  sophisti¬ 
cated  solvers  and  smoothing  techniques  and  hence  is  computa¬ 
tionally  cumbersome.  Fixed  grid  methods  on  the  contrary  [13-15], 
use  the  concept  of  similarity  transformation  to  convert  the  moving 
boundary  problem  to  fixed  coordinates.  One  such  method,  known 
as  Landau  transformation  [43],  is  used  in  the  solution  of  electro¬ 
chemical  model  for  lithium  ion  cells  and  is  adapted  in  this  work. 

In  context  of  LFP,  using  Landau  transformation  two  new  static 
position  variables  u  and  v  are  introduced,  one  for  each  phase. 
During  discharge  when  a  is  the  source  and  (3  is  product,  the  domain 
of  spherical  particle  Q  can  be  considered  to  be  constituted  of  QaU  Qp 
where  £aE  [0,  r“^p(t)]  and  Q^e[R,r^^(t)\.  Using  the  spatial  vari¬ 
able  u  which  is  given  as  u  =  r/r*^(t),  fixed  coordinate  equation 
for  diffusion  in  a  phase  or  the  core  is  obtained.  Similarly  using  ‘v’ 
defined  as  v  =  (r  -  the  corresponding 

fixed  domain  equation  for  diffusive  transport  in  the  shell  is  ob¬ 
tained.  The  transport  equations  during  charging  can  also  be 
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Fig.  3.  Voltage  vs  capacity  curve  prediction  for  different  current  densities  using  the 
discharge  model.  Also  comparison  of  present  approach  with  Srinivasan  and  Newman 
[15]  at  9  A  m-2  is  shown.  Result  from  Srinivasan  and  Newman  is  marked  by  legend  ‘S 
N’. 

obtained  using  similar  transformations.  For  continuous  charging 
after  discharging  (Fig.  2),  the  interface  formed  at  the  end  of 
discharge  will  define  the  boundary  condition  for  concentration  (Eq. 
(10)).  This  serves  as  the  initial  condition  for  charging,  when  new 
interface  (rj3^)  evolves  due  to  a  phase  formation.  In  order  to 
illustrate  this  point  a  schematic  of  charging  an  incompletely  dis¬ 
charged  spherical  particle  is  shown  in  Fig.  2.  The  figure  shows 
formation  of  multiple  cores  and  shells  but  it  is  important  to 
consider  that  only  the  outer  shell  of  a  phase  will  be  active,  and 
hence  only  will  govern  the  mass  transfer.  This  implies  tracking 
of  movement  of  a  single  interface.  On  the  contrary,  if  both  the 
interface  variables  (r?-^  and  r?^a)  were  used  to  track  the  dy¬ 
namics  of  phase  transition,  the  book  keeping  becomes  complicated, 
eventually  rendering  the  formalism  untenable  for  cyclic 
applications. 

Using  this  coordinate  transformation  the  equations  that  are 
specific  to  discharge  and  charge  that  evolve  out  of  the  general 
framework  (Eqs.  (6)-10)  are  listed  in  Table  2.  Although  the  number 
of  equations  increases  due  to  this  transformation,  the  formalism  is 
numerically  less  intensive  as  the  domain  does  not  change  with  time 
and  hence  can  be  adapted  to  simulate  charge-discharge  cycles.  The 
complete  numerical  solution  is  obtained  by  combining  the  equa¬ 
tions  given  in  Table  2  with  the  porous  electrode  model  (Eqs.  (1)— 
(5))  and  is  solved  using  a  Pseudo  2-D  [32-36]  approach. 

Numerical  implementation  in  a  multi-physics  environment  of 
COMSOL  [44]  is  obtained  by  suitably  modifying  the  available 
lithium  ion  cell  model  to  incorporate  multi-phase  coexistence.  This 
is  achieved  by  coupling  the  1-D  lithium  ion  cell  model  with  the  2-D 
solid  phase  diffusion  model  along  with  the  interface  movement 
equation.  Implicit  time  dependent  solver  with  BDF  (Backward  dif¬ 
ferentiation  formulas)  for  time  marching  is  chosen  with  an  absolute 
tolerance  of  1  x  10-3  and  relative  tolerance  of  1  x  10-2.  Initial 
numerical  studies  show  that  these  values  for  the  tolerances  are 
adequate  to  achieve  the  desired  numerical  accuracy  with  accept¬ 
able  computation  times.  To  obtain  the  numerical  solution  at  each 
time  step  a  solver  that  suitably  conditions  the  sparse  system  and 
optimally  allocates  the  memory  for  a  parallel  computing 
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Fig.  4.  Distribution  of  concentration  over  pseudo  2D  domain  of  positive  electrode  at 
t  =  3000  s  for  applied  discharge  current  of  3.6  A  m-2. 


architecture  (MUltifrontal  Massively  Parallel  sparse  direct  Solver, 
MUMPS)  is  used.  Further,  a  Newton-Raphson  based  technique  is 
used  to  obtain  solution  for  the  fully  coupled  multi-physics  PDEs. 

5.  Results  and  discussion 

The  generalized  moving  boundary  problem  developed  in  the 
earlier  section  is  used  to  simulate  the  response  of  LFP/C  cells  under 
various  scenarios.  In  this  section,  the  details  of  these  scenarios  are 
presented,  geometric  and  material  parameters  that  are  used  in  the 
simulations  are  listed  in  Table  3.  Srinivasan  and  Newman  [15]  and 
Thorat  25]  report  the  equilibrium  potential  for  anode  and  cathode 
respectively.  Details  of  these  potentials  are  given  in  Appendix  A, 
along  with  the  electrolyte  parameters.  Overview  of  the  case  studies 
discussed  in  this  section  is  as  follows: 

1 )  Proposed  numerical  implementation  is  validated  by  comparing 
model  predictions  with  the  shrinking  core  [15]  model  for 
discharge. 

2)  Generalized  moving  boundary  model  with  diffusion  in  a  phase 
is  used  to  simulate  the  discharge  and  charge  response  for 
various  current  rates.  Effect  of  diffusion  in  a  phase  is  illustrated 
using  the  spatial  and  temporal  concentration  profile  in  the 
positive  electrode. 

3)  Path  dependent  response  is  simulated  for  a  charge-discharge 
cycle.  This  simulation  tests  the  applicability  of  the  proposed 
approach  under  realistic  load  conditions,  in  addition  to 
analyzing  path  dependent  asymmetric  charge-discharge  of  LFP / 
C  cells. 

4)  Parametric  study  is  conducted  to  assess  the  applicability  of 
present  approach  as  a  potential  design  tool. 

Detailed  outcomes  of  the  individual  case  studies  are  discussed  in 
the  following  subsections. 

5.1.  Model  validation 

In  this  section  the  model  is  validated  by  comparing  with  results 
reported  by  Srinivasan  and  Newman  [15]  for  the  discharge  curve.  It 
can  be  seen  from  Fig.  3  that  the  model  compares  well  with  the 
published  results  for  an  applied  constant  current  of  9  A  m-2. 
Discharge  curve  is  simulated  for  a  range  of  current  densities  from 
3.6  A  m  2  to  36  A  m-2.  Model  outputs  for  different  current  densities 
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Fig.  5.  Distribution  of  interface  movement  over  pseudo  2D  domain  of  positive  elec¬ 
trode  at  t  =  3000  s  for  discharge  current  of  3.6  A  m-2. 


are  also  seen  in  Fig.  3.  The  results  show  expected  dependence  on 
increasing  current,  i.e.  drop  in  end  capacity  due  to  lower  utilization 
of  electrode  particles;  decrease  in  mid-plateau  potential  due  to 
higher  Ohmic  losses  and  increase  in  slope  of  the  mid-plateau  region 
[15]  due  to  higher  transport  losses. 

The  proposed  model  can  be  used  to  study  the  spatial  distribu¬ 
tion  and  evolution  of  Li  concentration  in  the  LFP  electrode.  The 
distribution  of  concentration  in  particles  along  the  thickness  of 
positive  electrode,  over  the  pseudo  2-D  domain  at  t  =  3000  s  is 
shown  in  the  surface  plot,  Fig.  4.  In  Fig.  4  and  in  all  subsequent 
surface  plots  of  the  pseudo  2D  domain,  the  x  axis  denotes  the 
thickness  of  the  electrode,  and  they  axis  the  radius  of  the  particle.  It 
may  be  noted  here  that  x  =  0  represents  the  separator  end  and 
x  =  Lp  represents  the  current  collector  end  of  electrode.  It  is 
observed  (Fig.  4)  that  though  diffusion  is  modeled  to  propagate  in 
radial  or  normal  direction  at  the  particle  scale,  the  concentration 
front  propagates  in  the  tangential  direction.  This  observation  is  in 
good  accord  with  experimental  [22,23  results  for  LFP.  Moving  one 
step  further  one  can  also  see  that  phase  interface  variable 
changes  from  almost  zero  at  the  cathode  current  collector  end  to  its 
initial  value  (=0.99  R)  at  the  separator  end  as  shown  in  Fig.  5.  Value 
of  close  to  zero  signifies  interface  reaching  the  center  of  the 
particle. 

It  is  interesting  to  see  in  Figs.  4,  5  that  there  is  higher  concen¬ 
tration  and  hence  a  higher  utilization  at  the  current  collector  end  in 


Fig.  7.  Comparison  of  evolution  of  interface  variable  f  for  constant  current  charge  and 
discharge  of  54  A  m-2. 


comparison  to  the  separator  end.  This  feature  is  specific  to  LFP 
electrodes,  as  an  opposite  trend  is  observed  [45]  for  other  elec¬ 
trodes.  The  reason  for  the  tangential  propagation  of  the  concen¬ 
tration  front  as  well  as  the  higher  utilization  at  the  collector  end  can 
be  ascribed  to  the  lower  electrical  conductivity  of  LFP  particle.  A 
higher  utilization  implies  a  higher  rate  of  surface  reaction.  An 
electrode  with  a  higher  electrical  conductivity  leads  to  a  design  case 
where  the  electrolyte  conductivity  is  limiting.  This  eventually  re¬ 
sults  in  a  higher  reaction  flux  at  the  separator  end  where  the  ions 
are  in  abundance.  In  the  case  of  a  porous  LFP  electrode  however,  the 
controlling  process  is  electrical  conductivity,  and  hence  the  trend  is 
reverse.  Implications  of  this  feature  in  optimal  cell  design  needs  to 
be  investigated.  It  is  clear  however,  that  any  increase  in  conduc¬ 
tivity  will  enhance  the  uniformity  in  exchange  current  density  and 
hence  the  utilization.  This  assertion  is  quantitatively  attested  by 
means  of  parametric  study  in  Section  5.4.3.  The  focus  of  this  section 
is  on  validating  a  discharge  model  for  cells  with  phase  change 
electrodes.  Further,  diffusion  in  a  phase  is  required  to  describe 
charge-discharge  response  as  discussed  in  Sections  3  and  4.  Re¬ 
sults  from  numerical  simulation  of  the  general  charge-discharge 
model  are  discussed  in  subsequent  sections. 


5.2.  The  charge-discharge  response  of  LFP/C  cells 

In  this  section,  the  results  of  the  generalized  moving  boundary 
model  to  represent  discharge  and  charge  behavior  are  presented.  It 


Fig.  6.  Voltage  vs  capacity  curve  for  charging  and  discharging  for  different  current  densities. 
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Fig.  8.  Distribution  of  concentration  over  pseudo  2D  domain  of  positive  electrode 
during  discharge. 


may  be  noted  that  while  charging  LFP  electrodes,  lithium  deficient 
a  phase  forms  the  shell  in  contrast  to  discharge.  In  the  present 
simulation,  for  diffusivities  of  a  and  (3  phase,  experimental  values 
reported  by  Safari  and  Delacourt  [26],  1.184  x  10-18  m2  s-1  and 
3.907  x  1CT19  m2  s  1  respectively  are  used.  Detailed  parametric 
dependence  on  diffusivity  is  presented  in  later  sections.  The  cell 
voltage  against  capacity  (defined  as  the  product  of  current  density 
and  the  time  of  charge/discharge)  for  charge  and  discharge  at 
various  current  densities  is  shown  in  Fig.  6.  Charge  voltage  profile 
exhibits  similar  trends  as  discharge:  increase  in  initial  voltage, 
decrease  in  end  capacity  and  increase  in  mid-plateau  slope  with 
increase  in  current  density.  Importantly,  on  comparing  voltage 
profiles  of  charge  and  discharge  (Fig.  6)  much  higher  utilization  is 
seen  during  charge  in  comparison  to  discharge  at  higher  current 
densities.  At  lower  current  densities  the  difference  of  utilization 
between  charge  and  discharge  is  lesser.  This  result  is  in  accord  with 
experimental  observations  by  Srinivasan  and  Newman  [29],  where 
it  is  reported  that  the  decrease  in  utilization  (or  end  capacity)  with 
increase  in  current  density  is  smaller  during  charge  as  compared  to 
discharge.  Hence,  it  is  evident  that  this  model  can  be  used  to 
represent  and  assess  the  asymmetry  in  the  charge  and  discharge 
response  of  LFP/C  cells. 

It  may  be  mentioned  that  although  discharge  and  charge  volt¬ 
ages  show  asymmetric  behavior  (Fig.  6),  these  results  are  obtained 
using  the  same  Ueqi  and  hence  the  asymmetry  is  of  transport  origin 
(Srinivasan  and  Newman  [29]).  A  further  cause  of  asymmetry  could 
be  due  to  difference  in  the  c*q  and  c%q  values  that  initiate  phase 
transformation  from  a  to  p  phase  and  vice-versa.  These  concen¬ 
trations  correspond  to  different  equilibrium  voltages  (Deq),  and 
over-potentials  (</>s  -  </>e  -  Heq )  resulting  in  different  rates  of  charge 
transfer  and  diffusion  processes  between  the  lithium  rich  and 
lithium  poor  phases.  It  can  be  seen  in  Fig.  6  that  at  low  rates  of 
charge  and  discharge  the  voltage  capacity  curve  is  nearly  sym¬ 
metric,  and  as  the  rate  of  discharge  or  charge  increases  the  asym¬ 
metry  increases.  In  order  to  explain  this  observation,  the  interface 
movement  variable  during  a  -*  p  (?  =  rj^/R)  and  p  -►  a 
(?  =  rf^a/R)  transitions,  normalized  with  radius  of  particle,  is 
plotted  at  a  high  current  of  54  A  m“2  in  Fig.  7.  As  seen  in  Fig.  7  the 
velocity  of  the  interface  is  different  between  charge  and  discharge, 
owing  to  difference  in  diffusivity  of  a  and  p  phases.  Importantly  the 
onset  of  phase  transformation  is  also  different  due  to  difference  in 
diffusivity  and  equilibrium  concentrations. 

Subsequently  the  effect  of  modeling  a  phase  diffusion  is  inves¬ 
tigated  by  assessing  the  concentration  distribution  in  the  positive 
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Fig.  9.  Distribution  of  concentration  over  pseudo  2D  domain  of  positive  electrode  at 
end  of  pulse  I. 


electrode  during  discharge  as  shown  in  the  surface  plot,  Fig.  8.  In 
the  snapshot  shown  in  Fig.  8,  the  concentrations  of  both  the  phases 
are  separated  by  a  front  that  eventually  propagates  from  the  correct 
collector  to  the  separator  end.  Thus  the  propagation  of  the  con¬ 
centration  front  (in  x  direction)  is  tangential  to  the  direction  of 
lithium  intercalation,  as  the  latter  occurs  in  the  y  (radial)  direction 
into  the  particles.  Thus,  this  surface  plot  shows  the  tangential  phase 
front  movement  as  well  as  the  separation  of  the  a  and  p  phases.  In  a 
recent  work  by  Fergusson  and  Bazant  [46],  this  phenomenon 
observed  in  LFP  particles  is  highlighted  well  and  the  same  is 
captured  through  non-equilibrium  thermodynamics  based  phase 
field  theory.  In  comparison,  the  present  approach  uses  a  compu¬ 
tationally  simpler  framework  of  porous  electrode  theory  to  arrive  at 
similar  conclusions.  It  is  to  be  mentioned  that  low  solid  diffusivity 
of  LFP  particle  promotes  the  tangential  phase  front  propagation 
during  charge-discharge  process. 

It  can  be  observed  from  Fig.  8  that  the  two  phase  coexistence 
causes  higher  concentration  gradients  due  to  the  difference  be¬ 
tween  the  equilibrium  concentrations  of  the  two  phases.  It  is  re¬ 
ported  [47]  that  high  concentration  gradients  due  to  phase 
separation  lead  to  build  up  of  high  stress  gradients.  These  stress 
gradients  induced  during  intercalation  process  result  in  increased 
crack  formation  [48]  in  electrodes.  The  stress  induced  cracks,  in 
turn,  lead  to  poor  electric  contact  and/or  capacity  fade  due  to  active 
Li  loss  by  fragmentation  or  isolation  of  electrode  particles  and 
eventually  to  cyclic  instability.  Thus  the  stress  in-homogeneities  in 
due  course  of  intercalation  cycles  can  be  detrimental  from  material 
mechanical  fatigue  point  of  view.  Such  scenarios  can  be  averted  by 
choosing  an  operating  window  in  predominantly  one  phase  region, 
or  tailoring  the  materials  for  a  narrower  equilibrium  concentration 
window  or  increasing  the  phase  transformation  sites.  The  proposed 
approach  can  be  used  to  model  such  aging  processes  on  incorpo¬ 
ration  of  suitable  coupled  stress  effects. 


5.3.  Path  dependent  behavior  during  charge— discharge  pulse 

In  the  experimental  study  by  Srinivasan  and  Newman  [29]  and 
Roscher  [49]  the  asymmetric  behavior  of  LFP  cathode  between 
charge  and  discharge  is  summarized.  It  is  concluded  that  at  the 
same  current  density  LFP  electrode  shows  higher  utilization  during 
charge.  This  in  turn  confirms  the  asymmetric  utilization  between 
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charge  and  discharge.  The  path  dependence  in  LFP  system,  how¬ 
ever,  arises  when  the  particles  are  not  fully  charged  or  discharged  - 
if  the  phase  transformation  in  any  one  direction  is  not  fully 
completed.  For  instance,  a  fully  charged  particle  (stage  I  of  Fig.  1) 
discharged  till  stage  III  and  charged  again  shows  a  different 
behavior  if  similar  state  of  charge  is  obtained  using  a  different  path. 
Such  path  dependence  is  characteristic  of  materials  that  exhibit 
phase  transition. 

In  this  section,  using  the  generalized  moving  boundary 
approach,  the  path  dependence  response  is  investigated.  It  is  to  be 
noted  that  the  differences  in  transport  and  kinetic  parameters 
during  charging  and  discharging  is  not  considered  in  this  study.  It  is 
interesting  to  investigate  if  the  path  dependent  response  is 
exhibited  even  without  this  information.  In  this  subsection  com¬ 
plete  voltage  profile  for  a  continuous  charge— discharge  cycle  is 
simulated  for  following  the  two  types  of  pulse  inputs: 

I.  Constant  current  1-C  discharge  for  9000  s  followed  by  rest  of 

5000  s  and  charge  at  1-C  for  9000  s  followed  by  rest. 

II.  Constant  current  7-C  discharge  for  500  s  followed  by  rest  of 

900  s  and  then  charge  at  7-C  for  600  s. 

Pulse  I  is  chosen  to  analyze  the  phase  separation  and  path 
dependence.  Pulse  II  shows  the  effects  of  rapid  charging  and  dis¬ 
charging.  Pulse  II  is  studied  to  demonstrate  efficacy  of  present 
approach  under  realistic  drive  cycles.  Concentration  distribution 
across  the  positive  electrode  under  pulse  I  during  discharge  is 
shown  in  Fig.  8,  and  at  the  end  of  the  pulse  in  Fig.  9.  The  path 
dependent  behavior  of  phase  change  electrode  as  discussed  in  Fig.  1 
is  quantitatively  obtained  in  Figs.  8,  9,  where  the  successive 
occurrence  and  de-occurrence  of  phase  separation  are  clearly  seen. 
Temporal  evolution  of  concentration  at  the  surface  of  positive 
electrode  particles  located  at  current  collector  end  is  shown  in 
Fig.  10.  The  discontinuous  evolution  of  concentration  during  a  -►  (3 
and  p  -►  a  phase  transformation  during  discharging  and  charging 
of  LFP  particle  is  demonstrated  in  Fig.  10.  The  initial  concentration 
evolution  during  the  a  ->  p  transition  is  shown  in  the  inset  of  the 
figure.  The  asymmetric  nature  of  utilization  in  LFP  cell  during 
continuous  discharge-charge  regime  of  pulse  I  is  also  seen  in 
Fig.  10.  Further,  some  amount  of  lithium  (denoted  as  AC  in  Fig.  10) 
remains  unused  after  one  charge-discharge  cycle.  This  unused 


Fig.  10.  Concentration  vs  time  profile  at  particle  surface  located  at  current  collector 
end,  showing  the  difference  at  the  end  and  start  of  pulse  I. 


Fig.  11.  Voltage  vs  time  during  discharge  and  charge  for  constant  current  of  25  A  m  2 
(pulse  II). 


lithium  accounts  for  the  observed  phenomenon  of  capacity  loss  in 
LFP  cells.  One  can  also  see  that  the  concentration  at  the  end  of 
charge  is  not  uniform  as  shown  in  Fig.  9.  Interestingly  recent  results 
on  X-ray  absorption  spectroscopy  (XAS)  by  Ouvrard  et  al.  50]  also 
show  such  heterogeneity  in  LFP  particles.  A  manifestation  of  such 
heterogeneity  is  a  reduction  in  amount  of  active  Li+  ion  during 
intercalation  cycle  of  the  cell  which  eventually  [51-54]  leads  to 
capacity  loss.  Such  a  capacity  loss,  as  it  does  not  involve  irreversible 
side  reactions  or  structural  damage,  can  be  recovered  by  long  time 
relaxations. 

In  order  to  further  demonstrate  the  efficacy  of  the  proposed 
approach  in  bringing  out  the  path  dependency,  simulation  under 
pulse  II  is  performed  with  different  extents  of  discharge  as  seen  in 
Fig.  11.  The  protocol  starts  from  a  fully  charged  cell,  the  varying 
extents  of  discharge  corresponds  to  a  difference  in  end  of  discharge 
SOC  of  0.2  between  the  pulses.  The  corresponding  voltage  during 
the  subsequent  charging  pulse  is  shown  in  Fig.  11,  and  a  detailed 
view  is  given  in  Fig.  12.  It  can  be  gleaned  from  the  results  that 
varying  extents  of  discharge  alters  the  subsequent  charging 
response.  It  needs  to  be  mentioned  that  this  behavior  is  observed 
irrespective  of  the  fact  that  the  end  of  discharges  correspond  to 


Time  (s) 

Fig.  12.  Voltage  vs  time  highlighting  the  effect  of  different  extents  of  discharge  on  the 
charging  profile. 
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Fig.  13.  Difference  in  surface  concentration  vs  time  at  current  collector  end  of  cathode 
for  different  discharge  times. 

SOCs  where  Ueq  is  almost  constant.  The  cause  of  this  feature,  thus, 
can  be  ascribed  to  difference  in  the  concentration  distribution  in 
the  positive  electrode  that  would  respond  differently  between 
discharge  and  charge.  A  reversal  of  current  loading  from  discharge 
to  charge  would  induce  a  reversal  of  phase  transformation, 
resulting  in  differences  in  the  concentration  profile.  Specifically,  on 
discharging  for  longer  times,  the  surface  concentration  of  cathode 
particles  at  the  end  of  discharge  are  farther  away  from  c%q  (Fig.  13). 
During  the  subsequent  charging,  these  values  equilibrate  back  to 
c^q  before  a  shell  of  a  phase  is  formed.  Thus,  the  charging  profile  in 
Fig.  12  varies  according  to  the  preceding  extend  of  discharge.  It  can 
be  also  seen  from  Fig.  12  that  the  time  to  achieve  the  constant 
voltage  region  (which  corresponds  to  the  two  phase  coexistence)  is 
longer,  if  the  extent  of  discharge  is  larger.  Thus  the  path  de¬ 
pendency  occurs  due  to  varying  relaxation  mechanisms  corre¬ 
sponding  to  the  operating  history.  Such  observations  due  to  path 
dependency  in  the  present  numerical  study  corroborate  well  with 
the  trends  observed  experimentally  [48]  on  LFP  based  cells.  The 
magnitude  of  voltage  difference  observed  for  different  instants  of 
discharging  in  the  present  simulations  is  in  accord  with  experi¬ 
mental  observations.  Further,  it  may  be  reiterated  here  that  the 
voltage  difference  is  obtained  with  the  same  values  of  diffusivity  for 
both  the  solid  phases;  the  impact  is  expected  to  become  more 
significant  with  incorporation  of  such  change.  Sensitivity  of 


diffusivity  on  the  path  dependence  will  be  studied  in  a  later  section. 
As  discussed  by  Ramana  [31]  and  Roscher  and  Sauer  [55],  another 
reason  of  path  dependency  in  LFP  cathode  is  due  to  inherent  hys¬ 
teresis  in  open  circuit  voltage  during  forward  and  reverse  trans¬ 
formation.  Incorporation  of  such  hysteretic  OCV  will  accentuate  the 
asymmetry  in  the  charge-discharge  response. 

5.4.  Parametric  study:  assessment  of  cell  design 

In  this  subsection  results  of  parametric  sensitivity  on  cell  per¬ 
formance  are  presented.  Aim  of  this  study  is  to  investigate  the 
correlation  of  the  extent  of  phase  transformation  with  macroscopic 
cell  response.  The  studies  in  this  section  suggest  that,  useful  design 
conclusion  can  be  derived  using  a  parametric  sensitivity  analysis.  In 
this  section,  effect  of  the  following  three  important  parameters  is 
studied: 

1.  Solid  phase  diffusivity 

2.  Radius  of  particle 

3.  Electronic  conductivity 

The  above  choice  of  parameters  is  inspired  from  various 
research  findings  that  have  identified  the  two  important  de¬ 
ficiencies  of  LFP  electrode  as  higher  transport  losses  and  low  rate 
capability.  To  overcome  these  deficiencies  the  natural  choice  is  to 
increase  the  diffusivity,  decrease  the  particle  radius  or  increase 
electronic  conductivity.  Analysis  in  this  section  aims  at  identifying 
the  critical  parameters  that  can  enhance  the  cell  performance. 

5.4.1.  Effect  of  solid  phase  diffusivity  on  discharge  behavior 

Diffusivity  values  of  both  the  a  &  [3  phase  are  set  to  be  equal  and 
is  varied  from  8  x  10  16  to  8  x  10-20  m2  s-1  and  the  discharge 
response  for  a  constant  applied  current  of  3.6  A  nrr2  is  shown  in 
Fig.  14.  The  corresponding  interface  movement  is  shown  in  Fig.  15. 
It  is  seen  that  higher  diffusivity  results  in  faster  interface  movement 
and  higher  capacity  (Figs.  14  and  15).  Eq.  (9)  shows  that  solid  phase 
diffusivity  is  directly  proportional  to  interface  velocity.  For  the 
lowest  value  of  diffusivity,  steeper  concentration  gradients  result  in 
quick  saturation  of  the  surface  concentration  of  the  cathode  par¬ 
ticles  which  in  turn  result  in  drop  in  the  cell  voltage  (Fig.  14).  These 
results  in  a  scenario  where  particles  of  the  positive  electrode  are 
not  utilized  completely  (Fig.  14).  This  study  suggests  that  the 
optimal  value  of  diffusivity  would  be  8  x  10-18  m2  s-1,  as  increasing 
diffusivity  beyond  this  value  enhances  the  voltage  and  capacity 
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Fig.  14.  Discharge  curve  for  different  values  of  solid  phase  diffusivity  for  constant  Fig.  15.  Evolution  of  normalized  interface  variable  f  for  different  values  of  solid  phase 
current  discharge  of  3.6  A  m-2.  diffusivity  at  3.6  A  m-2. 
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Fig.  16.  Voltage  vs  time  profile  during  discharge  of  pule  II  to  highlight  the  effect  of 
different  diffusivities  of  a  phase. 

response  marginally  (Fig.  14).  Interface  movement  (Fig.  15)  also 
shows  similar  trends  beyond  the  value  of  8  x  10-18  m2  s_1for 
diffusivity. 

5.4.2.  Effect  of  solid  phase  diffusivity  on  path  dependence  during 
continuous  discharge  and  charge  pulse 

This  case  study  highlights  the  importance  of  solid  phase  diffu¬ 
sivity  on  the  path  dependence  of  cell  response.  In  order  to  examine 
this  effect,  simulation  under  pulse  II  (described  in  Subsection  5.3)  is 
performed  for  different  diffusivities  of  a  phase  while  maintaining  a 
constant  diffusivity  of  the  p  phase.  The  cell  voltage  variation  with 
capacity  is  shown  in  Fig.  16  for  discharge  at  different  values  of  the 
diffusivity  ratio  y,  defined  as  y  =  Da/Dp.  During  discharge,  the 
system  is  predominantly  in  the  two  phase  coexistence  region,  and 
as  the  (3  phase  evolves,  the  mass  transport  is  governed  by  the 
interface  movement.  The  interface  movement  is  governed  by  the 
difference  in  diffusion  flux  of  the  two  phases  (Eq.  (9)),  and  is  slower 
for  higher  diffusivity  of  a  phase.  This  in  turn  results  in  higher 
concentration  gradients  in  the  shell  (p  phase)  and  higher  voltage 
drop  during  discharge.  Upon  successive  charging  (Fig.  17)  however, 
the  system  is  predominantly  in  the  single  phase  (p)  region  and  the 
expected  dependence  on  diffusivity  is  observed.  Thus  it  can  be  seen 
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Fig.  17.  Voltage  vs  time  profile  during  charge  of  pule  II  to  highlight  the  effect  of 
different  diffusivities  of  a  phase. 


that  multi-phase  coexistence  brings  about  the  effect  of  diffusivities 
of  both  the  solid  phases  on  the  cell  voltage  profile. 

5.4.3.  Effect  of  particle  radius 

One  of  the  suggested  ways  to  improve  the  rate  capability  or 
transport  limitation  in  LFP  electrodes  is  to  reduce  the  particle  size. 
The  discharge  curve  at  a  current  density  of  3.6  A  m-2  for  particle 
radius  decreasing  from  52  nm  to  13  nm  is  shown  in  Fig.  18.  It  can  be 
seen  from  Fig.  18  that  a  decrease  in  particle  radius  decreases  the 
diffusional  resistances,  and  for  the  same  applied  current,  results  in 
higher  discharge  voltage.  The  discharge  capacity,  however,  in¬ 
creases  marginally  on  reducing  the  particle  radius  to  13  nm,  as  seen 
in  Fig.  18.  The  interface  movement  at  these  particle  radii  (Fig.  19) 
also  does  not  show  any  significant  differences.  Interestingly  recent 
experimental  studies  by  Chueh  et  al.  [56]  confirms  the  findings  of 
present  simulation  results  wherein  it  is  reported  that  reducing 
particle  size  may  not  be  useful,  due  to  lower  diffusion  coefficient  in 
LFP.  Any  further  enhancement  could  only  be  achieved  by  inducing 
the  onset  of  the  phase  transformation  of  a  larger  number  of  parti¬ 
cles.  It  may  also  be  mentioned  that  particle  radius  reduction  does 
not  enhance  the  cell  performance  during  charging  either. 


Fig.  19.  Evolution  of  interface  variable  f  for  different  particle  radii  at  constant  current 
of  3.6  A  m-2. 
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Fig.  20.  Discharge  curve  for  different  values  of  electrical  conductivity  at  constant 
current  of  36  A  m-2. 

5.4.4.  Effect  of  electronic  conductivity 

As  mentioned  earlier,  one  of  the  important  limitations  of  LFP 
electrodes  is  the  low  electronic  conductivity.  This  causes  higher 
Ohmic  drops  in  the  electrodes  and  low  rate  capability  leading  to 
reduced  active  material  utilization.  The  present  case  study  aims 
at  elucidating  the  effect  of  electronic  conductivity  on  cell  voltage 
and  discharge  capacity.  To  highlight  this  effect,  discharge  at  a 
higher  constant  current  density  of  36  A  itT2  is  performed  for 
three  different  values  of  conductivity  varying  in  multiples  of  1, 10 
and  100  from  the  base  value  of  5  x  10  3  S  itT1.  As  seen  in  Fig.  20 
the  corresponding  discharge  curve  shows  enhanced  performance 
with  higher  discharge  voltage  and  increased  utilization  with  in¬ 
crease  in  conductivity.  Increase  in  voltage  is  due  to  reduced 
Ohmic  drop  and  increase  in  utilization  (or  end  capacity)  is  due  to 
lower  transport  losses.  It  is  also  observed  from  the  concentration 
distribution  at  the  electrode  level  (figure  not  shown)  that 
reduction  in  Ohmic  drop  increases  the  reaction  flux  at  the  par¬ 
ticle  surface.  This  results  in  the  concentration  front  emanating 
from  both  the  collector  and  the  separator  end  of  the  electrode 
leading  to  higher  utilization.  It  is  to  be  noted  that  further 
increasing  the  conductivity  from  0.05  S  m-1  to  0.5  S  m_1  en¬ 
hances  the  performance  marginally  (Fig.  20)  and  hence  an 
optimal  choice  of  conductivity  is  0.05  S  m-1  for  better  cell  design. 
It  can  be  concluded  in  this  section  that  out  of  the  three  param¬ 
eters  studied  conductivity  has  the  dominant  effect  in  improving 
the  performance  of  LFP  based  cells. 


6.  Summary  and  conclusion 

In  this  work  a  generalized  moving  boundary  framework  for 
lithium  ion  cells  with  electrodes  that  exhibit  phase  transition  is 
proposed.  This  framework  is  used  to  develop  electrochemical 
model  based  on  pseudo  2-D  approach  of  porous  electrode  theory. 
Using  suitable  coordinate  transformation  the  moving  domain 
problem  is  translated  in  to  an  equivalent  fixed  coordinate  system. 
Suitable  path  dependent  algorithm  is  devised  to  model  the  elec¬ 
trochemical  response  under  continuous  charge  and  discharge 
cycle.  Numerical  simulations  show  efficacy  of  the  present 
approach  in  suitably  capturing  the  underlying  electrochemical 
processes  in  LFP/C  cells  namely,  tangential  front  propagation, 
phase  separation,  charge-discharge  asymmetry  and  path 
dependent  response.  The  concentration  front  propagates 
tangentially  in  due  course  of  discharge  owing  to  lower  diffusivity 
and  electrical  conductivity  of  LFP  electrodes.  Generalized  core 


shell  model  ascribes  sharp  interface  formation  between  the  two 
solid  phases  of  LFP  which  adequately  captures  the  phase  separa¬ 
tion.  Phase  separation  in  LFP  is  a  very  important  physical  phe¬ 
nomenon  as  it  causes  steep  concentration  gradients.  These 
gradients  result  in  high  stress  generation  and  are  tantamount  to 
mechanical  fatigue.  Charge-discharge  study  at  various  currents 
predicts  asymmetry  in  cell  utilization,  and  especially  at  higher 
rates  it  is  found  that  utilization  is  higher  during  charge  in  com¬ 
parison  to  discharge.  This  prediction  of  lower  utilization  at  higher 
current  density  during  discharge  as  compared  to  charge  corrob¬ 
orates  well  with  the  experimental  results.  Comparison  of  interface 
movement  during  charge  and  discharge  confirms  the  difference  in 
utilization  and  ascribes  this  asymmetry  to  be  of  transport  origin. 
Another  important  feature  of  the  proposed  generalized  frame¬ 
work  lies  in  its  ability  to  model  the  path  dependent  response. 
Robustness  of  present  numerical  implementation  is  shown  by 
predicting  the  response  for  both  high  and  low  rate  charge- 
discharge  pulse  protocol.  Also  it  is  shown  that  the  charging 
voltage  profile  depends  on  the  extent  of  preceding  discharge, 
resulting  in  the  path  dependent  response  of  phase  change  LFP 
material.  Parametric  sensitivity  studies  suggest  that  among  the 
important  properties  of  LFP  (diffusivity,  particle  radius  and  elec¬ 
tronic  conductivity)  increase  in  conductivity  enhances  cell  per¬ 
formance  significantly.  A  higher  electrical  conductivity  enhances 
the  phase  transformation  resulting  in  to  improved  cell 
performance. 

The  generalized  moving  boundary  framework  proposed  here 
can  be  coupled  with  energy  balance  equations  to  develop  a  non- 
isothermal  model  to  facilitate  thermal  management  studies. 
Assessment  of  carbon  coating  or  solid-solution  doping  and  other 
enhancement  features  can  be  incorporated.  Furthermore  the  model 
can  be  coupled  to  mechanical  momentum  balance  equations  to 
predict  stress  gradient  due  in  a  phase  change  electrode  material  to 
develop  fatigue  models  to  study  battery  aging.  Some  of  these  ex¬ 
tensions  are  in  progress. 

Appendix  A 


This  appendix  enlists  the  equilibrium  potential  for  LFP  cathode 
which  is  given  in  Ref.  [25]  as: 


U  =  2.56742  +  57.69(  1  —  tan  h  (  100 - +  2.9163927 

Cmax 


0.442953  tan-1  (  -  65.41928 


Qnax 


0.097237  tan-1  (  -160.9058 


64.89741 


154.59 


Cmax 

Similarly  the  equilibrium  potential  for  graphitic  anode  is 
adopted  from  Ref.  [16]  and  is  given  in  terms  of  electrode  state  of 
charge  (x)  as: 

U(x)  =  0.124+ 1.5  exp(-70x)) -0.0351  tan  h  (-  °;°286) 

\  0.083  J 


-  0.0045  tanh 
-0.0147  tan  h 


x  -  0.9 
04A9 
x  -  0.5 
0.034 


-  0.035  tan  h 

-  0.102 tanh 


x  -  0.99 
005^ 
x- 0.194 
0.142 


-  0.022  tan  h  ( *  -  0.01 1  tan  h  ( *Q 


0.0155  tan  h 


0.0164  ) 
x- 0.1 05 
0.029 


The  electrolyte  properties  -  diffusivity  and  conductivity  15], 
used  in  the  present  study  are  given  respectively  as: 
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D  =  5.34  x  l(T10exp(  -  1.565  x  lO^CejmV1  s  source 

V  /  p  product 

Ke  =0.0911  +  1.901  x  10_3ce  -  1.052  x  10“6ce 

+  0.1554  x  lO^CeSm-1  Subscript 

app  applied 

where  ce  is  the  electrolyte  concentration  mol  m-3.  disch  discharge 


List  of  symbols 


as  Active  surface  area  per  electrode  unit  volume  for  electron 

transfer  reactions,  m-1 
Brug  Bruggeman  coefficient 

c  lithium  concentration,  mol  m-3 

cs  lithium  concentration  at  the  surface  of  electrode,  mol  rrr3 

Co  initial  lithium  concentration,  mol  m  3 

ce  Li+  concentration  in  the  electrolyte,  mol  m-3 

cseq  equilibrium  lithium  concentration  for  species  ‘s’,  mol  m  3 
cmax  maximum  lithium  concentration  in  LFP  particles, 
mol  rrr  3 

C  current  needed  to  completely  discharge  the  electrode  in 

an  hour 

D/  diffusion  coefficient  of  Li+  in  phase  i,  m2  s_1 

De  diffusion  coefficient  of  the  electrolyte,  m2  s-1 

f±  mean  molar  salt  activity  coefficient 
io  exchange  current  density,  A  rrr2 

I  discharge  current,  A 

io  Exchange  current  density,  A  rrr  2 

in  intercalation  current  density,  A  m-2 

k  kinetic  rate  constant 

Q.  discharge  capacity,  A  h 

rza^  position  of  the  phase  boundary  when  transforming  from 
a  to  p  phase,  m 

rf^a  position  of  the  phase  boundary  when  transforming  from 
P  to  a  phase,  m 

Rg  gas  constant,  8.3145  J  mol-1  K-1 

R  radius  of  LFP  particles,  m 

5  geometric  area  of  the  electrode,  m2 

SOC  state  of  charge 

t  time,  s 

t®  transference  number  of  the  electrolyte 

T  temperature,  K 

U  equilibrium  potential  of  the  electrode,  V 

be  Boundary  condition 

PDE  Partial  Differential  Equation 


Greek  letters 

a  Lithium  poor  phase  of  LiFePCU 

p  Lithium  rich  phase  of  LiFePC^ 

oca>  occ  Transfer  coefficient 

4>e,  4>s  liquid  or  solid  phase  potential,  V 

5P>  5S  electrode  or  separator  thickness,  m 

re  Volume  fraction  of  electrolyte 

Ke  Conductivity  of  the  electrolyte,  S  m-1. 

<7 s  Conductivity  of  the  solid  phase,  S  rrr1. 

Q  Physical  domain  of  a  or  p  phase  in  the  particle 

vs^p  Velocity  of  interface  for  transforming  from  source  to 
product 


Superscript 

ref  reference  state  at  298.15  I< 

eff  effective  property 
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